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Abstract 



The apparent stability of population oscillations in ecological systems is a long-standing 
puzzle. A generic solution for this problem is suggested here. The stabilizing mechanism 
involves the combined effect of spatial migration, amplitude-dependent frequency, and 
noise: small differences between spatial patches, induced by the noise, are amplified by 
the frequency gradient. Migration among desynchronized regions then stabilizes the oscil- 
lations in the vicinity of the homogenous manifold. A simple model of diffusively coupled 
oscillators allows the derivation of quantitative results, like the dependence of the desyn- 
chronization upon diffusion strength and frequency differences. For an unstable system, 
a noise-induced transition is demonstrated, from extinction for small noise to stability if 
the noise exceeds some threshold. The coupled oscillator model is shown to reproduce all 
previously suggested stabilizing mechanisms. Accordingly, we suggest using this model 
as a standard tool for identification and classification of population oscillations on spatial 
domains. 
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I. INTRODUCTION 



The apparent stability of prey-predator systems is an age-old puzzle. While natural 
host-parasitoid and prey-predator systems persist for millions of years without extinction, 
simple considerations imply that this situation should be inherently unstable. Essentially, 
when a predator consumes a prey, it clearly increases its fitness and its chance to breed 
and produce more offspring. Thus, the predation rate should behave autocatalyticaly, 
leading to the extinction of the prey species and, consequently, the predator population. 
Even before the information and quantitative era, ancient day naturalists were able to 
recognize this puzzle. Among them, Herodotus and Cicero perceived the persistence of 
prey species in the face of adversity as a manifestation of divine power and the creator's 
design (Cuddington, 2001). Herodotus' explanation for the phenomenon was based on 
multiplication rates. In his words: "... for timid animals which are a prey to others are 
all made to produce young abundantly, so that the species may not be entirely eaten up 
and lost; while savage and noxious creatures are made very unfruitful." This explanation, 
clearly, can not stand on its own. Even if Herodotus' statements about the reproduction 
habits of certain species were true (he claims that the hare is able to superfetete, i.e., 
to conceive while still pregnant; while the lioness may use her womb only once, since 
the cub ruptures it with his claws during pregnancy), balancing two exponentially di- 
vergent processes (multiplication and predation) is achievable only via fine-tuning of the 
parameters. Moreover, that fine-tuning must be robust against various kinds of noise. 
Indeed, as a rule, natural prey-predator systems are always subject to noise: either some 
sort of environmental stochasticity (e.g., years with more/less rain, weather fluctuations, 
heterogenous habitat) or, even under fixed conditions, some noise is introduced into the 
system due to demographic stochasticity, related to the discrete nature of the individ- 
ual agents and the stochastic character of birth-death-predation processes (Lande, 2003; 
Foley 1997). 

In modern times, the course of explanations has been changed, becoming more deter- 
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ministic and quantitative. The basic idea relies on the fact that when the prey population 
decays there is not enough food for the predators and their population also diminishes, 
so both populations oscillate around their mean values. As pointed out by Nicholson 
(1933), these oscillations are an intrinsic property of interacting populations. If, for 
example, the density of a host is above its steady value, it will be reduced by the para- 
site. However, when the host reaches its steady density, the density of parasites is above 
its steady value. "Consequently, there are more than sufficient parasites to destroy the 
surplus hosts, so the host density is still further reduced in the following generation ... 
Clearly, then, the densities of the interacting animals should oscillate around their steady 
value." This paradigm was formulated mathematically, using deterministic continuous- 
time partial differential equations (PDE's), by Lotka and Volterra (LV model) (Lotka, 
1920; Volterra 1931; Murray 1993). An analogous model with a discrete time step was 
introduced for a parasitoid-host system by Nicholson and Bailey (NB) (Nicholson and 
Bailey, 1935). Both models [and their variants (May, 1978; Hassell and Varley, 1969; 
Crowley, 1981; Rosenzweig and MacArthur, 1963)] allow for a coexistence fixed point, i.e., 
a steady state that admits finite population of the predator and the prey, and predict 
population oscillations around this coexistence steady state. The prediction of popula- 
tion oscillations is one of the major achievements of the deterministic models. Its most 
well known example is century-long documentation of oscillations in the Canadian lynx - 
Snowshoe Hare inhabitants of Northern America (Elton, 1924; Stenseth et al., 1998). 

However, beyond this success, Lotka, Volterra and Nicholson recognized that their 
scheme fails to solve the basic stability problem, as it cannot account for the sustain- 
ability of these oscillations (Nicholson, 1933; Cuddington, 2001). If, under the influence 
of demographic or environmental stochasticity, oscillation amplitude increases, it should 
eventually reach an extinction point for one of the species. But this is exactly the situa- 
tion for systems described by the Lotka- Volterra model: as oscillations of any amplitude 
are allowed, any type of noise will eventually increase their amplitude and the species 
will go extinct (see section 2 below). For the Nicholson-Bailey map, the situation is even 
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worse: the coexistence fixed point is not only marginally stable but unstable. Thus, even 
without a continuous effect of noise, any small perturbation should increase exponentially 
and drive the system to extinction. Trying to maintain the coexistence state in a NB 
host-parasite model is like trying to balance a cone on its tip. 

Facing the inability of the LV and NB models to admit oscillation sustainability, one 
may feel a temptation to abandon them altogether in favor of another formulation that 
supports an attractive fixed point or (if oscillations are required) a limit-cycle/strange- 
attractor. While it is reasonable to assume that some natural populations are described 
by such models, experimental and field studies show that, quite generically, this instability 
is a main characteristic of reality. In fact, for small systems and well- mixed populations, 
experiments show oscillation growth and extinction in a wide variety of interacting species 
(Gause , 1934; Pimentel et al., 1963; Huffaker, 1958). Spatially-extended systems, on 
the other hand, seem to support sustained oscillations, as emphasized by field studies, 
experiments, (Kerr et al., 2002; Kerr et al., 2006; Luckinbill, 1974; Holyoak and Lawler, 
1996) and numerics (Wilson et al., 1993; Bettelheim et al., 2000; Washenberger et al., 
2006; Kerr et al., 2002; Kerr et al., 2006). 

These findings give an appeal to Nicholson's (1933) old proposal about migration in- 
duced stabilization, i.e., that desynchronization between weakly coupled spatial patches, 
together with the effect of migration, stabilize the global populations. To get the flavor of 
the mechanism, let us imagine a metapopulation on two patches, where within each patch 
the population oscillations are governed by, say, Lotka-Volterra dynamics with constant 
diffusion (i.e., constant per capita migration rate) between these two patches. Clearly, as 
emphasized in Figured], if the oscillations on these two patches desynchronized, e.g., if 
one of the patches is in a state of large populations while the other is, at the same time, 
in a diluted state, migration between patches pushes the whole system inward toward 
the coexistence fixed point yielding sustained oscillations. However, one should bear in 
mind that migration is a double edged sword, as it tends to avoid population gradients 
and leads to synchronization. The acid test for Nicholson's proposal is thus as follows: is 
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FIG. 1: Population oscillation on two spatial patches coupled by migration. If both patches 
desynchronize, one may find one of them (A) in the dense population state and the other one 
(B) in the dilute phase. Diffusion tends to decrease population gradients, hence the whole system 
flows towards the coexistence fixed point, represented in the lower panel by an asterisk. 

the diffusion among patches weak enough to allow noise-induced desynchronization, but 
at the same time strong enough to stabilize desynchronized patches? If this is to be the 
case, the desynchronization-diffusion stabilization may work. 

Unfortunately, examination of this idea in many studies, summarized in a recent review 
article (Briggs and Hoopes, 2004), yields negative results. Generically, diffusion stabilizes 
the homogenous manifold and different spatial patches get synchronized, leading back 
to the well-mixed unstable dynamics (Crowley, 1981; Allen, 1975; Reeve, 1990). At the 
end of the day, thus, we are back where Herodotus left us more than 2500 years ago: we 
cannot explain the sustainability of prey-predator systems, despite their presence since 
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the beginning of life on earth. The understanding of oscillation persistence is not only an 
intellectual challenge; stabilization of such oscillations is considered to be a major factor 
affecting species conservation and ecological balance (Earn et al., 2000; Blasius et al., 
1999). 

The goal of this paper is to provide a generic solution to this puzzle, a solution that 
is independent of external assumptions, such as space-time fluctuations or heterogenous 
migration patterns. We will show that the basic ingredient that leads to desynchroniza- 
tion is the dependence of oscillation frequency on their amplitude (Amplitude Dependent 
Frequency, ADF), and will survey the consequences and the implementation of this idea 
for different systems and types of noise. The basic quantitative and qualitative tool pre- 
sented here is a simple toy model for coupled oscillators. Apart from its usefulness in the 
analysis of the ADF mechanism, this toy model works for all types of stabilizing mech- 
anisms presented so far, and serves as a prediction and classification tool. As such, we 
suggest our coupled oscillators system to be used as the standard -modeling technique for 
metapopulation oscillations. 

The present work is organized as follows: in the next section, the stability problem will 
be introduced for a single patch system, together with a short review of the experimental 
and numerical results that support the instability for well-mixed populations. In the 
following section, Nicholson's proposal - desynchronization induced stability on spatial 
domains - is endorsed, based again on numerics and experimental results. The fourth 
section is the main part of the paper, where a toy model is presented and analyzed 
and the role of amplitude-dependent oscillations is clarified. Next we give a short critical 
review of previous solutions to the stability problem, and show how to incorporate them all 
within the framework of diffusively-coupled oscillators. We end up with some concluding 
remarks, trying to sketch a classification scheme for observed population dynamics in 
experimental and numerical situations. 
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II. INSTABILITY AND EXTINCTION OF WELL MIXED POPULATIONS 



We begin the demonstration of the instability problem looking at the Lotka-Volterra 
predator-prey system, the paradigmatic model for oscillations in population dynamics 
(Lotka, 1920; Volterra, 1926; Murray, 1993). The model describes the temporal evolution 
of two interacting populations: a prey (b) population that grows with a constant birth rate 
a in the absence of a predator (the energy resources consumed by the prey are assumed 
to be inexhaustible), while the predator population (a) decays (with death rate /x) in the 
absence of a prey. Upon encounter, the predator may consume the prey with a certain 
probability. Following a consumption event, the predator population grows and the prey 
population decreases. For a well-mixed population, the corresponding partial differential 
equations are: 

— = — iia + X\ab (1) 

L/C 

^ = ab — X 2 ab 
at 

where Ai and A2 are the relative increase (decrease) of the predator (prey) populations 
due to the interaction between species, correspondingly. 

The system admits two unstable fixed points: an absorbing state a = b = and 
the state a = 0, b = 00. There is one coexistence, marginally stable fixed point at 
a = a/X 2 , b = /i/Ai. Local stability analysis yields the eigenvalues ±i^/JIa for the sta- 
bility matrix. Moreover, even beyond the linear regime there is neither convergence nor 
repulsion. Using logarithmic variables z = ln(a), q = ln(b) eqs. ([I]) take the canonical 
form z = dH/dq, q = —dH/dz, where the conserved quantity H (in the ab representa- 
tion) is: 

H = Xib + X 2 a — /i ln(a) — a ln{b). (2) 

The phase space is thus segregated into a collection of nested one-dimensional trajectories, 
where each one is characterized by a different value of H, as illustrated in Figure [2j Given 
a line connecting the fixed point to one of the "walls" (e.g., the dashed line in the phase 
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FIG. 2: The Lotka-Volterra phase space (left panel) admits a marginally stable fixed point 
surrounded by close trajectories (three of these are plotted). Each trajectory corresponds to a 
single H defined in Eq. ([2]), where H increases monotonically along the (dashed) line connecting 
the center with the a = wall, as shown in the lower right panel. In the upper right panel, the 
period of a cycle T is plotted against H, and is shown to increase almost linearly from its initial 
value T = 2ir / \/JIa close to the center. 

space portrait, Figure [2]), if is a monotonic function on that line, taking its minimum 
Hmin at the marginally stable fixed point (center) and diverges on the wall. It turns out 
that all the important features related to the instability and the synchronization depend 
only on the topology of the phase space, and the actual values of the growth, death and 
predation parameters are not important. Thus, without loss of generality, we employ 
hereon (unless otherwise stated) the symmetric parameters fi = a = \i = \ 2 = l. The 
corresponding phase space, along with the dependence of H on the distance from the 
center and a plot of the oscillation period vs. H, are represented in Figure [2j 

Given the integrability of that system, the effect of noise is quite trivial: if a and b ran- 
domly fluctuate in time, the system wanders between trajectories, thus performing some 
sort of random walk in H with "repelling boundary conditions" at H min and "absorbing 
boundary conditions" on the walls (as negative densities are meaningless, the "death" of 
the system is declared when the trajectory hits the zero population state for one of the 
species, i.e., when at least one of the species becomes extinct). Thus, under the influence 
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FIG. 3: The survival probability Q(t) is plotted versus time for a single patch noisy LV system. 
Eqs. ([T]) (with the symmetric parameters) were integrated numerically (Euler integration with 
time step 0.001), where the initial conditions are at the fixed point a = b = 1. At each time step, 
a small random number r](t)At was added to each population density, where rj(t) £ [—A, A]. 
A typical phase space trajectory, for A = 0.5, is shown in the inset. The system "dies" when 
the trajectory hits the walls a = or b = 0. Using 300 different noise histories, the survival 
probability is shown here for A = 0.5 (full line), A = 0.3 (dotted line) and A = 0.25 (dashed 
line). Clearly, the survival probability decays exponentially at long times, Q(t) ~ expi—t/r), as 
expected for a random walk with absorbing boundary conditions, 1/r scales with A 2 . 

of noise, thus, the problem is reduced to that of a random walk with an absorbing trap [a 
first passage problem (Redner, 2001)]. While the first passage problem is characterized by 
power law decay, here the topology of the orbits forces the trajectory back to the vicinity 
of the absorbing walls, and hence the decay is exponential. An example of a phase space 
trajectory for a single patch noisy LV system is shown in the inset of figure [3X and the 
average growth of the oscillation amplitude is clearly seen. The survival probability Q(t) 
(the probability that a trajectory does not hit the absorbing walls within time t) is shown 
for different noise amplitudes in figure [3j 

At this point it is instructive to clarify the types of noise to be used throughout 
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this paper. The first and the most simple is the additive noise, i.e., random addition or 
subtraction of population density at certain time steps. In the real world, this corresponds 
to migration of small parts of the population in and out of the considered system. As the 
implementation of an additive noise in Langevin type equations is relatively simple, we 
use it as the "noise of choice" for the present work. 

Indeed, there are two common sources for time-dependent randomness in ecological 
systems. One source is the environmental stochasticity, i.e., fluctuations in the environ- 
mental conditions that may change the reaction and interaction parameters. The growth 
rate of the prey, for example, should be taken not as a constant a but as a fluctuating 
quantity around some average, (a) + 5a(t), where (• • •) stands for average over time. 
Another source of noise is the intrinsic demographic stochasticity. While Eqs. (JT]), for 
example, describe continuous populations sizes a(t) and b(t) that may take any value, the 
actual number of individuals in a community should be an integer (Lande et al., 2003). 
Indeed, Partial differential equations like ([I]) only describe the dynamic of the average 
population and are derived from the exact master equation as a "mean field" approxi- 
mation, neglecting higher correlations and approximating, e.g., (a(t)b(t)) ~ (a(t)) (b(t)) . 
The actual dynamics of the system, described by the master equation, is not determinis- 
tic but stochastic, and this implies that noise of order is applied to any population 
of iV individuals. While the relative amplitude of that noise decays like 1 / yN for large 
populations, it is still of importance if there is no attractive manifold for the deterministic 
dynamics (Kessler and Shnerb 2006), such as in the LV or NB models. The effect of demo- 
graphic stochasticity on a single path LV system is very similar to that of additive noise, 
as indicated by the simulation results obtained using Gillespie's event driven algorithm 
(Gillespie, 1977). 

In the current work, additive noise is used in order to model multiplicative noise such 
as demographic stochasticity. This procedure may be justified using a self-consistency 
argument: we want to present a mechanism that stabilizes the population oscillations, 
such that the number of individuals in, say, the predator population, does not deviate 
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strongly from its average value. If this is the case, the 0(yN) noise amplitude does 
not change so much along the orbit and the system "feels" constant (additive) noise. 
The noise amplitude of environmental stochasticity is proportional to the population size 
[the term 5a (t) presented above multiplied by the prey population in ([I])], thus it is also 
approximated by an additive noise with a constant amplitude if the oscillations are not 
too large. 

While the LV system is marginally stable and is driven to extinction by the noise, the 
situation for the host-parasitoid model of Nicholson and Bailey is even worse. In their 
model, the host population H and the parasite population P satisfy the map (Nicholson 
and Bailey, 1935): 

H t+1 = aH t e- xp \ 
P t+1 = cH t (l - e- XPt ) (3) 

where a > 1 is the growth factor of the host in the absence of a parasite. The system 
admits a coexistence fixed point at: 

H '"" 



c\(a — 1) 

and small deviations from that fixed point grow exponentially with time, yielding a prompt 
extinction of one of the interacting species (see Figure H]). Here the effect of noise is 
negligible: any infinitesimal deviation from the fixed point renders the system extinction- 
prone. 



III. SPATIAL STRUCTURE AND SUSTAINED OSCILLATIONS 

As concluded in the last section, a well-mixed population that obeys Lotka-Volterra 
or Nicholson-Bailey dynamics is extinction prone: oscillation amplitude grows in time 
until one of the species goes extinct. The only difference between the two models is the 
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FIG. 4: Nicholson-Bailey dynamics. The map ([3]) is iterated 90 times where the initial deviation 
from the fixed point is 0.001. The parameters are a = 1.2, A = 0.3, c = 6.1. 

way this extinction is approached: while an NB system is unstable and oscillations grow 
in time exponentially even without noise, the LV model admits neutral oscillations and 
approaches extinction only due to the effect of stochasticity. 

One may thus suggest that something is wrong in the modeling of interspecific in- 
teraction. While LV and NB are clearly the simplest descriptions of the corresponding 
systems, maybe in nature some more complicated processes take place, and more accu- 
rate mathematical models should use equations that support an attractive manifold, i.e., a 
phase space region that "attracts" the trajectories. In particular, one may write a model 
that supports a stable fixed point, limit cycle (Kolmogorov, 1936; May, 1972) or strange 
attractor. In the first case, after a short period of oscillations the system settles into a 
stable state of constant populations; in the second, it will converge to fixed amplitude 
oscillations; and in the third scenario (strange attractor), it will wander chaotically be- 
tween a set of points confined to some phase space region. Any of these solutions is robust 
against small noise, as the dynamic is dissipative, small fluctuations flow back into the 
attractive manifold. 

However, a series of experimental results shows that, indeed, for small sample sizes, 
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prey-predator or host-parasitoid systems actually flow to extinction. From the seminal 
work of Gause on the Didinium - Paramecium system (Gause, 1934), through Huffaker's 
orange experiment (six spotted mite - Typhlodromus ) (Huffaker, 1958), and Pimintel's 
wasp-fly setup (Pimentel et al., 1963); all of these experiments demonstrated that popu- 
lation dynamics result in extinction of one of the species. 

Of particular interest is a new series of experiments dealing with the cyclic, rock-paper- 
scissors dynamic (Kerr et al., 2002; Kirkup and Riley, 2004; Kerr et al., 2006) . In the 
first set of experiments, colicin-producing strains (C) of E. Coli were interacting with 
colicin-sensitive (S) and colicin-resistent strains (R); while C kills S, S outcompetes R 
and R wins against C. It turns out that for well-mixed systems, both in vitro (Kerr et al., 
2002) and in vivo (Kirkup and Riley, 2004), fluctuations drive two of the three species 
to extinction via the same mechanism of "random walk" among naturally stable orbits 
(Reichenbach et al., 2006). In the second experiment (Kerr et. al., 2006) E. Coli and 
its parasitoid (phage) were subjected to controlled migration. Again, in the "well-mixed" 
regime, when all the biological materials were mixed each 12 hours, extinction took place 
after a very short time. 

On the other hand, the spatial structure and, in particular, spatial migration have 
been shown to be of crucial importance in the sustainability of interacting populations. 
The famous field studies about natural prey-predator systems [like the Canadian Lynx - 
Snowshoe Hare (Elton, 1924; Stenseth et. al., 1998) data from Hudson Bay company, the 
Moose- Wolf coexistence in Isle Royal (Wilmers, et al., 2006), and the biological control of 
the Prickly pear cactus by the moth cactoblastis cactorum in eastern Australia (Freeman, 
1992)] indicate that, in spatially-extended systems, the population oscillations are either 
negligible or finite, ensuring persistence. Luckinbill's (Luckinbill, 1974) experiment, using 
again the Didinium - Paramecium system, shows that the system's lifetime (time until 
extinction of one of the species) grows almost exponentially with its linear size. Similar 
observations were made by Holyoak and Lawler (1996), and, as already mentioned, by the 
series of experiments dealing with rock-paper-scissors dynamics (Kerr et. al., 2002; 2006) 
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, where slow mixing yields finite fluctuations and promotes survival. There are also many 
numerical experiments on stochastic models (like the individual-based models, subject to 
demographic stochasticity) that show stable oscillations for large samples(Wilson, et al., 
1993; Bettelheim et al, 2000; Washenberger et al., 2006 ;Kerr et. al, 2002; 2006). 

All in all, much evidence supports the idea that interacting species dynamics are in- 
herently unstable on a single patch, and that the stability observed in nature is related to 
the fact that natural systems are spatially-extended. This idea may be implemented in an 
extinction-recolonization scenario (Taylor, 1990), but in general this mechanism works for 
any system that supports desynchronization. As pointed out by Nicholson (1933), if the 
system desynchronizes, the migration between patches will induce an inward flow toward 
the coexistence fixed point and stabilize the amplitude of oscillations. The main puzzle, 
thus, is reduced to the identification of the conditions under which desynchronization 
takes place. It turns out, however, that generically the diffusion tends to lock the system 
to its synchronized state, stabilizing the homogenous manifold such that the whole system 
acts like a single patch, rendering unbounded oscillations (Briggs and Hoopes , 2004). 

The condition for desynchronization in diffusively coupled patches have been examined 
in many studies and the main results, summarized in a recent review article (Briggs and 
Hoopes , 2004), are as follows: 

• For any network of N patches, if the migration between patches is symmetric or 
almost symmetric (i.e., the diffusion of the prey and the predator are, more or less, 
the same), there is no diffusion-induced instability, and the homogenous manifold 
is stable (Crowley, 1981; Allen, 1975; Reeve, 1990). Thus, the effect of migration 
alone does not cure the instability problem. 

• The system may become desynchronized when in the presence of spatial heterogene- 
ity, e.g., where the reaction parameters vary on different spatial patches (Murdoch 
and Oaten, 1975; Murdoch et al., 1992; Hassell and May, 1988). In that case, the 
intrinsic dynamic at any localized patch takes place on different time scales for 
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the same concentrations, so diffusion fails to synchronize different patches. This 
mechanism may be generalized to include not only "quenched" heterogeneity but 
also environmental stochasticity, i.e., where the reaction parameters are subject to 
spatio-temporal fluctuations (Crowley, 1981 ; Reeve, 1990; Taylor, 1998) . 

• Spatial heterogeneity may be introduced into the model in a more sophisticated way: 
one may define a homogenous dynamics in all patches with an underlying nonlin- 
ear process that supports spontaneous symmetry breaking and pattern-formation. 
These patterns, in turn, will control the reaction parameters to yield desynchroniza- 
tion. 

• Diffusion-induced instability may occur if the migration rate of the predator is much 
smaller than that of the prey, particulary if the prey migration rate is zero [Jansen, 
1995; Abta et al. (2006 b)]. 

While spatial heterogeneity, environmental stochasticity, and different migration rates 
for different species are indeed characteristics of many natural systems, these mechanisms 
are not generic and it seems that they cannot explain the experimental results that show 
persistence in "clean" extended systems, like (Luckinbill, 1974) and (Kerr et al.; 2002; 
2006). Accordingly, it seems that the combined effect of noise and diffusion is a necessary 
precondition for population stabilization. However, up until now, the qualitative nature of 
the "missing" underlying mechanism has remained obscure, and no theoretical framework 
that allows for quantitative prediction has been presented. 

This theoretical gap has been addressed recently (Abta et al., 2006 a) using a toy model 
for coupled oscillators presented in the next section. The main new ingredient empha- 
sized by the proposed model is the dependence of frequency on the oscillation amplitude, 
reflected by the gradient of the angular velocity along the radius u'(r). The instabil- 
ity induces desynchronization iff the small, noise-induced differences between patches are 
amplified by the frequency gradient such that the "desynchronization parameter" (<p 2 ) 
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acquires a finite value, leading to "restoring force" toward the origin of the homogenous 
manifold. 

Moreover, it turns out that our toy model may also reproduce the other mechanisms 
suggested in the past, and yield simple analytical results together with qualitative and 
quantitative predictions. In the fifth section, these generalizations are presented. Thus, 
our toy model allows for a comprehensive typification of different mechanisms for sustained 
oscillation on diffusively coupled patches, and yields simple criteria for classification of 
experimental results and identification of the stabilizing mechanism. 



IV. TWO PATCH SYSTEM 



To consider the effect of spatial structure, let us begin with the simplest case, i.e., two 
patch system connected diffusively as in Figure [TJ One should bear in mind, though, that 
any finite system that admits an absorbing state and is subject to additive noise must 
eventually reach extinction. Under the influence of noise, any stable state becomes only 
metastable, as some rare configuration of the noise will push one of the species to extinc- 
tion. The appearance of sustained oscillations, i.e., of an attractive manifold, manifests 
itself in finite system only in its lifetime before it hits the absorbing state. This lifetime, r, 
is known in the literature as the "intrinsic mean time to extinction" (Grimm and Wissel, 
2004), and if a metastable state appears, it increases exponentially with the "strength" of 
attraction (the Lyapunov exponent) of that manifold. If the stability of a system increases 
with the number of patches, one may expect the lifetime of the system to diverge with 
the number of patches, yielding real stability for an infinite system. The best example of 
this phenomenon - a metastable state in a small system becoming stable as its size grows 
- is the logistic growth with demographic stochasticity, where the transition from extinc- 
tion to proliferation belongs to the directed percolation equivalence class (Grassberger, 
1982; Janssen, 1981). A similar transition, belonging to the same universality class, was 
observed recently for coupled limit cycles with demographic stochasticity (Mobilia, et al, 
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2006). Thus we are not trying to show that the two patch system is really stable; expo- 
nential growth in the lifetime of the system is a clear indication of the appearance of an 
attractive manifold that yields real stability in the infinite size limit. 
The simplest example is the LV system on two patches: 



dt 
da 2 

~dt 
dbi 

~dt 



—fiai + Aiai&i + D a (a 2 — ai) 
-fia 2 + Xia 2 b 2 + D a (a 1 - a 2 ) (5) 
abi - X^bi + D b (b 2 - 61) 



Ob? 

— = crb 2 - \ 2 a 2 b 2 + D b (b x - b 2 ). 

The time evolution of that system, with additive noise, equal diffusivities, D a = D b = 
D, and symmetric reaction rates is obtained through Euler integration. In the limit D = 
the patches are unconnected; thus, starting from the homogenous fixed point a\ = a 2 = 
b\ = b 2 = 1, the single patch situation still holds and the system hits the absorbing 
walls after a characteristic, noise-dependent, time. In the opposite limit, D = 00, the 
system sticks to the invariant manifold and acts like a single patch (with modified noise 
and interaction parameters), again performing a random walk in the invariant manifold. 
However, between these two extremes, there is a region where the combined effect of 
diffusion and noise stabilizes a finite region within the invariant manifold. For any noisy 
two-patch system Q(t), histograms (like those shown for a single patch in Figure EJ) were 
used to extract the typical decay time ["intrinsic mean time to extinction" (Grimm and 
Wissel, 2004)] t(D,A) by fitting its tail to exponential decay exp(—t/r). In Figure [3, 
t(D, A) is plotted against D for different noise amplitudes, and is shown to increase 
(faster than exponentially) with D (1/D) as D approaches zero (infinity). Evidently, for 
intermediate diffusivities, an attractive manifold appears in phase space, with a Lyapunov 
exponent that grows (faster than linearly) with the diffusion constant. 
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FIG. 5: The typical persistence time as a function of the diffusion rate for different levels of 
noise. The values of r were gathered from survival probability plots (like those in Figure [3]) and 
are displayed here for the two-patch system. One sees that the value of r grows very rapidly 
(even faster than exponentially) with the migration rate for small diffusion values, and decays 
with D for large diffusivities. Data is shown for different noise intensities A = 0.3 (triangles), 
0.5 (squares) and 1.0 (circles). 

A. The Coupled oscillator model and LV system 

Although both the Lotka-Volterra and Nicholson-Bailey models are oversimplified with 
respect to the complexity of interspecies interactions in the natural world, they are still too 
complicated to allow for an understanding of the stabilization mechanism when migration 
and noise interfere. The main obstacle (as will become clear later) is that the radial 
velocity along a trajectory depends not only on the amplitude of oscillations (related 
to the conserved quantity H), but also on the location along a trajectory, i.e., on the 
azimuthal angle between the point along the trajectory and the fixed point 9 (see Figure 
E]). In order to clarify the origin of the stable cycles, let us introduce a toy model that 
imitates the main features of the real systems. Although that model does not allow for 
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predator 

FIG. 6: The angular velocity along some orbits of the Lotka-Volterra dynamic. Fast regions are 
marked in red, slow regions in blue. Clearly, the dynamic is slowest when the populations of 
both species are diluted, and fastest along the dense region in the upper-right "shoulder." Note 
that the velocity gradient along an orbit increases with H. 

an absorbing state, it captures the basic mechanism for stabilization of spatially-extended 
systems in the presence of noise. 

The toy model deals with the phase space behavior of diffusively-coupled oscillators, 
where the angular frequency depends on the radius of oscillations. With additive noise, 
the Langevin equations take the form, 

dx\ 



u){x x ,yi)yi + D x {x 2 - x x ) +T]i(t) 
uj(x 2 ,y 2 )y 2 + D t (xi - x 2 ) + r] 2 (t) 



dt 
dx 2 

~dt 

—1 = + D 2 (y 2 - y x ) + r} 3 (t) 

dy 2 
dt 



(6) 



-lu(x 2 , y 2 )x 2 + D 2 (y x - y 2 ) + r) 4 (t), 



where all the 7/-s are taken from the same distribution. If the angular frequency is loca- 
tion independent, u(x,y) = u , the problem is reduced to coupled harmonic oscillators, a 
diagonalizable linear problem that admits two purely imaginary eigenvalues in the invari- 
ant, homogenous manifold. With noise, the random walk on that manifold is independent 
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of the motion in the fast manifold and the radius of oscillation diverges with the square 
root of time. 

Now let us define the oscillation radius for each patch, = \Jxf -\- yf for i — 1,2, 
and assume that the angular frequency depends only on that radius and is ^-independent 
[9i = arctg(yi/xi)}. With that, the total phase $ = 61+62 decouples and the 3-dimensional 
phase space motion is dictated by the equations (we take D 1 = D 2 = D and define 
= Q x - e 2l R = n + r 2 , r = r x - r 2 ): 

R = -2Dsin 2 (<p/2)R + f) R (7) 
r = -2Dcos 2 {<p/2)r + fj r (8) 

As before, all the fj-s are taken from the same distribution. Note that <fi represents the 
phase desynchronization between patches while r is the amplitude desynchronization. 

Eqs. ([7]) clarify the role of phase desynchronization as the stabilizing mechanism. The 
dynamics in the homogenous (R) manifold look very much like that of an overdamped 
harmonic oscillator in noisy environment, 

z = -kz + T](t), (10) 

a well-known system [Ornstein-Uhlenbeck process, (Gardiner, 2004)] that admits the 
steady state Boltzman distribution P(z) ~ exp(—kz 2 /A 2 ). However, if the phases of 
these two patches synchronize and the expectation value of 2 vanishes, so does the 
"spring constant" of the oscillator. Without that restoring force, the motion on the R 
manifold is a simple random walk, so the oscillation amplitude grows indefinitely. Phase 
{(f)) desynchronization, thus, is the crucial condition for stabilization. This feature is 
stressed in the inset of Figure [7J where the flow lines of the deterministic dynamics in the 
R — (j) plane are sketched: the line = is marginally stable, but any deviation leads to 
inward flow. 
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Close to the invariant manifold, when and r/R are much smaller than one, the 
amplitude desynchronization r is solvable. Neglecting corrections of order 2 , 



and the r 2 typical fluctuation around zero is of order A 2 /D. However, the amplitude 
desynchronization factor r does not appear in ([7]), and the stabilization is determined only 
by the phase. Accordingly, the system supports an attractive manifold iff the amplitude 
desynchronization yields phase desynchronization. 

Another piece of information may be gathered from the harmonic limit, u(r) = uq. 
Here there should be no phase synchronization, as we already diagonalized the linear 
equation and find no restoring term in the homogenous plane. Looking at Eq. ([9]) with 
uj{ri) —uo{r2) = one concludes, thus, that the rightmost (noise) term in (J9]) is irrelevant. 
The dependence of the frequency on the amplitude (i.e., the dependence of u on r) should 
be the factor that allows the translation of amplitude desynchronization into a phase 
desynchronization. Intuitively, when two patches with different oscillation amplitudes 
move with different angular velocities, this immediately yields phase differences. 

Without the explicit noise term, Eq. (Q may be written as, 



where the approximation is valid close to the invariant manifold. Again we face an over- 
damped harmonic oscillator, where now the source of noise is the r fluctuations (obeying 
the Boltzman statistics). With that, 



f = —2Dr + 7], 




which is again an equation for an overdamped harmonic oscillator, so 




(12) 




(13) 




(14) 



may be plugged into ((Tj) to yield: 




(15) 
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FIG. 7: Histograms showing the probability to be at distance R from the origin as a function 
of R, for two coupled noisy oscillators, where u = 1 + ar with D = 0.01 with various values 
of noise strength A and angular velocity gradient a. As expected, the phase space confinement 
is proportional to a, from a = 1 (triangles) to a = 0.5 (circles) to a = 0.1 (solid line), all 
for the same level of noise A = 0.1. On the other hand, as predicted by the linear analysis 
close to the invariant manifold, the confinement is noise-independent, and the three solid lines 
corresponding to different levels of noise (A = 0.1, 0.5, 1) with the same a = 0.1 almost coincide. 
The inset shows the flow lines on the r = plane. The invariant manifold = is stable, but 
there is a "maximally desynchronized" unstable orbit converging to the center at <j) = ir. If the 
expectation value of 4> 2 deviates from zero, there is an effective restoring force toward the center, 
and the noise-induced random walk on the = manifold is bounded. 

Now the "restoring force" along the R coordinate (the invariant manifold) is finite, and, 



This radius of stable oscillations diverges as D — > oo, as expected. The small D instability 
(decoupled patches) manifest itself in the divergence of r 2 as D — > 0. Surprisingly, since 
both the restoring force and the noise in the invariant manifold are proportional to A 2 , 
the expected R distribution has to be noise-independent at that limit, as demonstrated 
in Figures [71 and H~2l A more careful examination of the predictions about the expectation 



(i? 2 ) 



D 2 



(16) 
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FIG. 8: ((f) 2 ) (triangles) (r 2 ) (squares) and (R 2 ) (circles) as a function of a, for two coupled 
oscillators where uj = 1 + ar, D = 0.01 and A = 0.1. Clearly the value of (r 2 ) is almost 
independent of the frequency gradient a. (R 2 ) is fitted to a -1,99 , in agreement with Eq. ( I16p . 
{(ft 2 ) scales like a 1,74 , close to the prediction of (JHJ). 




FIG. 9: ((j) 2 } (triangles) (r 2 ) (squares) and (R 2 ) (circles) as a function of the migration rate 
D, for two coupled oscillators where u> = I + ar, a = 0.4 and A = 0.1. As predicted, (r 2 ) is 
inversely proportional to D (best fit to L> -0 - 98 ). (R 2 ) is fitted to D~ L88 , [predicted exponent is 
(-2)] and ((f) 2 ) scales like D~ 2 52 (predicted value -3). 
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D=0.01 a=0.1 




FIG. 10: (eft 2 ) (triangles) (r 2 ) (squares) and (R 2 ) (circles) as a function of the noise Delta, 
for two coupled oscillators where uj = 1 + ar, a = 0.1 and D = 0.01. As predicted, (r 2 ) is 
proportional to A 2 (best fit to A -1,98 ), ((f) 2 ) is fitted to A 1 ' 47 , (predicted exponent is 2) and 
(R 2 ) clearly saturates at large noise levels. 

value of the three parameters (Eqs. [TBI HH an d the Boltzman distribution for amplitude 
desynchronization) is presented in Figs. El El and dOj The small deviations from the 
predicted exponents are related to the failure of the approximations used far from the 
invariant manifold. Decreasing the noise level, these exponents approach their predicted 
values. Figure [TT1 shows the D dependence of the lifetime for the coupled oscillator model; 
its agreement with figure [5] is evident. 

Another direct manifestation of the coupled oscillator predictions in the LV system 
is presented in Figure [12j Again, for "intermediate" migration (e.g., D = 0.01), the 
average distance from the origin saturates, while the chance to find the system at large H 
becomes exponentially small, as illustrated in Figure [121 In agreement with the results of 
the toy model, the flow toward the center is correlated with the phase desynchronization, 
leading to stabilization of oscillations at finite H. As predicted, while the width of the <j) 2 
distribution depends strongly on the noise amplitude, the oscillation amplitude is almost 
noise independent. 
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FIG. 11: The typical persistence time as a function of the diffusion rate for different levels of 
noise for the coupled oscillator model. The "extinction" was declared when the amplitude of 
oscillations for one of the patches exceeded unity. The values of r were gathered from survival 
probability plots (like those in Figure [3|). One clearly observes the similarity with Figure 
Data is shown for different noise intensities A = 0.3 (triangles), 0.4 (squares) and 0.5 (circles). 
Note that while the expectation value of R 2 is noise independent, the noise still controls the 
probability to rare fluctuations that leads to extinction. 

B. The Coupled oscillators model and NB instability 

In order to imitate the behavior of the NB system with its unstable fixed point, Eqs. 
([6]) should be modified to include this new ingredient, manifested here by the terms 
proportional to k, 

kxi + lu(xi, y x )yi + D x (x 2 - x x ) + 771 (t) 

KX 2 + UJ(X 2 , y 2 )V2 + D X (X X - X 2 ) + 7] 2 (t) 

Ky x - uj(x u yi)x! + D 2 (y 2 - y x ) + rj 3 (t) (17) 




dx\ 

~dt 

dx 2 

~dt 
dyi 
dt 
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FIG. 12: The average AH at an elementary time step (0.001 of a unit time) as a function of 
the angle (j) between the patches. While a simple phase space random walk yields, on average, 
positive AH, this property is shown here to hold only for small 4>. At larger angles, the diffusion 
between patches forces the system toward the center and the average AH becomes negative. 
Results are shown for A = 0.1 (full line) and A = 1 (dashed line). The inset shows the 
probability distribution function for H at these two noise levels. 

= ny 2 -uj(x 2 ,y 2 )x 2 + D 2 (y 1 - y 2 ) +rj 4 (t). 

The system is still invariant with respect to global rotation, thus it reduces to the three 
dimensional phase space: 

R = ( K K-2Dsin 2 ( ( f)/2)) R + fi R (18) 
r = (k -2D 'cos 2 '((f)/ '2)) r + fj r (19) 

* = "2 D {^Pj sin<t> + W (r 2 ) - ^ + (| - |) . (20) 

Close to the invariant manifold, thus, 

P(r) ~ exp[-(2D - k))t 2 / A 2 } (21) 

and 



J (r) 2 A 2 
,0(2.0 - k) 2 ' 
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FIG. 13: Extinction times for the unstable, coupled oscillator cartoon of the Nicholson-Bailey 
model. Eqs. (|17p were integrated numerically (using Euler integration with At = 0.002) for 
K = 0.0001 and and uj = + for different noise amplitudes A and diffusion constant D. The 
extinction time is plotted for four different noise level against the diffusion constant, and the 
two transitions are implicit. The lifetime of the system for large noise (A = 0.1,0.2, triangles) 
diverges beyond our computational abilities for D = 0.01. 

Consequently: 

2 A 2 A 2 

( R f ~ £(02) _ - = uj\r) 2 A 2 /{2D - n) 2 - n ^ 

The system becomes unstable when either r 2 or R 2 diverge. The first criterion for 

stability comes from the amplitude synchronization parameter, 2D > k, so the diffusion 

should increase above some threshold value in order to prevent desynchronized extinction 

where the system acts as if made of two disconnected patches. As in the LV case, if 

the migration rate is too large (i.e., if k becomes larger than [u (r)] 2 A 2 /(2D — k) 2 ~ 

u' 2 A 2 /D 2 ), the system synchronizes and the deterministic flow leads to synchronized 

extinction. Note that close to the low D transition the extinction rate grows with the noise, 

while close to the second transition, increase of noise amplitude A yields lower extinction 

rates, emphasizing the fact that the stability is noise-induced. This feature is clearly seen 

in Figure [TU where for small noise the oscillation amplitude grows exponentially in time 
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FIG. 14: Noise-induced transition for the coupled oscillator cartoon of the Nicholson-Bailey 
model. Eqs. (|17p were integrated numerically (using Euler integration with At = 0.002) with 
k = 0.0001, D = 0.01 and uj = ujq + r/2 for different noise amplitudes A. The total distance 
R (averaged over 100 runs) is presented, in logarithmic scale, against time measured in units of 
ujq. Small noises are followed by exponential growth of the oscillation amplitude, as suggested 
by the deterministic part of (|17p . The larger the noise, the slope of this diverging line becomes 
smaller. If the noise is large enough, R saturates at a finite value, as seen more clearly when the 
scale is not logarithmic (inset). 

while for large noise the oscillation stabilizes. 

Along this section, the various numerical implementation of the toy assumes, for sim- 
plicity, that the oscillation frequency depends linearly on the slope. If, on the other hand, 
uj' is a nonuniform function of the amplitude, oscillations will grow sublinearly (for the 
LV model) or exponentially (in NB case) with time, until they reach a phase space region 
where duj/dr is large enough. For an LV-like system it leads to the appearance of a "soft" 
limit cycle - exponential convergence from the outside, diffusion inside. 

It should be noted that the stability mechanism presented here is applicable not only 
for the classic predator-prey systems but for any biological system that supports, locally, 
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neutrally stable or unstable oscillations. In particular it may stabilize a " locally" unstable 
ecology of interspecific competition for common resource. Such a system may be described 
mathematically by the generalized Lotka-Volterra equations (Chesson, 2000); even if the 
theory predicts extinction, species diversity may be maintained on spatial domains. An- 
other possible application is the stabilization of "eigenperiods" - oscillations due to the 
interaction between an individual and its internal resources or by maternal effect - sug- 
gested recently as a predation free mechanism for population oscillations (Ginzburg and 
Colyvan, 2004). 



V. OTHER MECHANISMS 



This section is devoted to a discussion of previously suggested stabilization mechanisms 
in the framework of the coupled oscillator toy model. We show that all the significant 
mechanisms surveyed by (Briggs and Hoopes , 2004) become very clear and allow analyt- 
ical understanding when presented using the coupled oscillator system. 



A. Spatial heterogeneity 

To implement the spatial heterogeneity (SH) (Murdoch and Oaten , 1975; Murdoch 
et al, 1992; Hassell and May, 1988) in our framework, consider a two patch system with 
different (radius independent) frequencies, iO\ and u 2 : 

dx\ 



dt 

dx 2 

~dt 



viVi + D 1 (x 2 - xi) + ?7i(t) 
uj 2 y 2 + D 1 (x 1 - x 2 ) + 772(0 



^ = - UlXl + D 2 (y 2 - Vl ) + m (t) (24) 

-J^ = -u 2 x 2 + D 2 (yi -y 2 ) +m( t )- 

As these equations are linear, one may again diagonalize the deterministic part of the 
evolution matrix to obtain the exact probability distribution function. However, as this 

30 



problem is still invariant under global rotation and is independent of Q\ + 6 2 we have 
chosen to use the more intuitive polar representation. After the standard coordinates 
transformation, one finds: 



R = -2Dsin 2 R + fj R (25) 
r = -2Dcos 2 (J^j r + rj r (26) 

<P=-2D + [u 2 - <*] + (| - |) • (27) 

Here Au = u 2 — oj\ is a constant, so close to the invariant manifold, r decouples 
from 0. Consequently, while P(r) is still given byO satisfies an equation for a forced 
overdamped pendulum and sin(<j)) = Au/2D. The "spring constant" in the invariant 
manifold R is thus noise-independent and 

<* 2 > - (S- (28) 

The radius of oscillation in a heterogenous system increases with noise: the system be- 
comes more extinction prone as the noise grows larger. In the corresponding "Nicholson- 
Bailey" toy model the phase transition happens at k = D/(Auj) 2 , and the location of that 
transition is noise-independent. 



B. Environmental stochasticity 

The above framework may also be used to consider the stabilizing effect of spatio- 
temporal environmental stochastisity (STES) (Crowley, 1981; Reeve, 1990; Taylor, 1998). 
If on both patches the radial velocity takes the form u + Q(t), where i is the patch index 
and ( is a white noise that satisfies (C(^)C(^O) = Y5(t — t'), the <fi variable obeys an 
equation for an overdamped noisy pendulum and the resulting motion on the invariant 
manifold satisfies: 

R = (k-T)R + Vr , (29) 
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so (R 2 ) ~ A 2 /(T — k). In that case the lifetime of the system is diffusion-independent, 
growing with the environmental stochasticity and decaying with the noise. 



C. Jansen instability 

About ten years ago, Jansen (Jansen, 1995; Jansen and de Roos, 2000; Jansen and 
Sigmund, 1998) put forward the idea of linearly unstable orbits of the Lotka-Volterra dy- 
namics, i.e., orbits in the homogenous manifold for which the highest absolute value of an 
eigenvalue of the Floquet operator is larger than 1. This may happen only if the migration 
properties of the prey and the predator are different. In the case of equal diffusivities, 
the migration term factors out from the Floquet operator and the stability properties of 
orbits lying in the invariant manifold are the same as their matching trajectories on a 
single patch (Abarbanel, 1995). However, Jansen pointed out that if one sets = in 
Eqs. El some orbits may become unstable. In that case one may use the fact that the 
total H 

H T = Hi + H 2 = a\ + a 2 + &i + b 2 - Zn(aia 2 &i&2) , (30) 

with the deterministic dynamics (j5J) is a monotonously decreasing quantity in the non 
negative population regime: 

= _ Da a^t\ _ Dt nj^m < (31) 

at \ aia 2 ) \ b x b 2 j 

Accordingly, if an orbit on the invariant manifold becomes unstable, the flow will be 
inward and the population oscillations stabilize. 
With the transformation, 

A = ^l B= b -1±^ (32) 
2 2 v ; 

a_ 1 -a 2 ^b^-h 
2 2 ' 

one realizes the homogenous AB manifold and the 5 — 6 coordinates that measures the 
deviation from that manifold (the heterogeneity of the population). In these coordinates 
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the system satisfies, 

BA 

— = -fiA + X^B + XtSe (33) 

— = aB- X 2 AB + X 2 86 
Bt 

= -u5 + \ 1 AB + \ l B5-2D a 8 

at 
BO 

— = oB- X 2 AO - X 2 B5 - 2D b 9. 
Bt 

Linearizing around the homogenous manifold, the AB dynamic is equivalent to that of a 
single patch, 

A = —jj,A + X t AB (34) 
B = aB - X 2 AB 



and the 5 — 6 linearized dynamic is 
/ 




fi + X 2 A - 2D a X 2 

-X 2 A a-X 2 B- 2D b 



(35) 



V 



One may thus calculate the eigenvalues of the Floquet operator for one period along an 
orbit of (!34|) . The resulting stability diagram, first obtained by (Jansen, 1995), is shown 
in Figure [151 

We have discussed Jansen's stabilization mechanism in a different publication (Abta 
and Shnerb, 2006b). It turns out that the underlying mechanism has to do with the 
dependence of the angular velocity along the orbit on the azimuthal angle (see Figure [6]), 
and we can mimic that behavior using our toy model with u{6). Specifically, the coupled 
oscillator model with, 

Bx\ 



Bt 
Bx 2 

~BT 



u{0x)yx + D x (x 2 - xi) 
u{9 2 )y 2 + D x (x x - x 2 ) 



(36) 
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FIG. 15: Stability diagram for phase space orbits (ordered by their conserved quantity H) for 
different values of predator diffusion D a , where Db=0. 



I = -wftjxi + D y (y 2 - y x ) 



dy_ 
dt 

—1 = - u ( y 2 ) X2 J r D y (y l -y 2 ), 



where D x = D, D y = and 



uj = uj + ui cos{6 - -J, 



(37) 



leads to the same type of instability. 

To prove that this model actually yields Jansen's instability, we have used (i G 1,2) 



r- = X 2 , + yj 



6i = arctan( — 



(38) 



. _ (xx + yy) g _ (xy- yx) 



and 



r = r 2 — r 1 

<p = e 2 -e 1 



R = r 2 + ri 
$ = 9 2 + X . 



(39) 
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The flow in the invariant manifold satisfies, 

R = 

6 = w(0i)+w(0 2 ), 



(40) 



while the linearized equations for the desynchronization amplitude r and the desynchro- 
nization angle <fi satisfy: 



20^ sin<E> 



2u/($/2) -4D x cos 2 ($/2) 

-D a . J Rsin($) -4D x sin 2 ($/2) 




(41) 



Using the Floquet operator technique to analyze the stability of an orbit by integrating 
ETj) along a close trajectory of (j4"U|) . one finds the stability map presented in Figure 
[TBI where the parameter if of Figure [TH] is now replaced by u\, which measures the 
"eccentricity" of the angular velocity along a circular path. Here, two unstable regions 
appear, for large and small D x , but the qualitative picture is almost the same. 



VI. CONCLUDING REMARKS: TOWARDS CLASSIFICATION OF SUS- 
TAINED OSCILLATIONS. 

Up until now, four different mechanisms that induce stability of population oscilla- 
tions for metapopulations have been presented. The generic mechanism relies on the 
dependence of the oscillation frequency on their amplitude (ADF, Amplitude Dependent 
Frequency). The other three are spatial heterogeneity (SH), spatio-temporal environmen- 
tal stochasticity (STES), and Jansen's mechanism based on the dependence of the angular 
velocity on the azimuthal angle. 

In natural systems, in the laboratory, and in simulations, one may encounter population 
oscillations in the presence of a few of the above mentioned factors, simultaneously. The 
task of the researcher is to make a distinction between them and to identify the relevant, 
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FIG. 16: Stability diagram in the u\ — D x plane for the Floquet operator (same as Figure [T5|) 
for the coupled oscillator system described by Eqs. (|36p with D y = 0. The two unstable regions 
correspond to different signs of the Floquet unstable eigenvalue, as explained in the text. 

or the most relevant, stabilizing mechanism. To do that, the following comments may be 
useful: 

• Some of the parameters involved, such as the migration rate per capita, the strength 
of the environmental stochasticity, growth and competition parameters and so on, 
may be known from the literature or may be measured independently. If possible, 
an estimate of the stable oscillation amplitude has to be made and a comparison 
between possible mechanisms will reveal the dominant one, which is (in the long 
run) the one that stabilizes the smallest oscillation amplitude. 

• If the dominant mechanism is spatial heterogeneity, the phase between cycles on 
neighboring patches is kept constant in time. Thus, it is quite easy to recognize an 
SH-induced stability. In all other cases there is no preference among patches and <fi 
is distributed randomly around zero. 

• Jansen's RVA instability depends very much on the difference between species dif- 
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fusivities. It is impossible to observe such an instability in host-parasite systems, 
where the parasite migration depends on host movement. 

• One way to distinguish between STES and ADF is to separately manipulate the 
additive noise (taking small portions of the population in and out of the system 
at random) and environmental stochasticity. An increase in the amplitude A of 
the additive noise leads to larger oscillations for STES and smaller/equal oscillation 
amplitudes for ADF; amplified stochasticity will diminish oscillations if the basic 
mechanism is STES, and is neutral in ADF-controlled systems. Another test is the 
response of the system to changes in the migration rate, where STES synchronized 
oscillation amplitudes are not affected by changes in D. 

• In both SH and STES mechanisms, the amplitude desynchronization decouples from 
the phase desynchronization, while in ADF they are intimately connected. Measure- 
ments of the correlations between the phase and the amplitude desynchronization 
will immediately reveal the relevant mechanism. 

As pointed out by Jansen and Sigmund (Jansen and Sigmund, 1998), "all models of 
ecological communities are approximations: it is pointless to burden them with too many 
contingencies and details. On the other hand, they would be of little help if they were not 
robust against the kind of perturbation and shocks to which a real ecosystem is ceaselessly 
exposed." The current accuracy of data on population oscillations and inter-specific in- 
teraction, both from field studies and from experiments, is, as far as we know, far beyond 
the level needed for an exact comparison with theoretical predictions about the oscillation 
phase portrait, like those predicted by Lotka-Volterra and other models. Given that, the 
main insights from the available data are, first, the mere existence of these oscillations, 
and second, the identification of the underlying mechanism that limits the amplitude of 
these oscillations in noisy environments. As emphasized by the recent experiments of 
(Kerr et al., 2002 ; Kerr et al., 2006) , one may observe persistent oscillations or extinc- 
tion, but it is hard to compare the exact population dynamic with the predictions of the 
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theory. Accordingly, the analysis of data on population cycles may be preformed, as we 
have shown here, completely within the framework of the simple coupled oscillation model 
that allows all the suggested limiting processes within a transparent and general modeling 
scheme. 
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